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Integrateurs symplectiques en geometrie sous-riemannienne: 

le probleme de Martinet 



Resume : On compare les performances d'integrateurs symplectiques et non symplectiques pour le calcul de 

geodesiques normales et de points conjugues dans un example sous-riemannien, le probleme de Martinet. On 
etudie le probleme d'abord avec une metrique plate, puis avec une perturbation a un parametre conduisant a 
des geodesiques non integrables. De cette etude, on deduit que proche des directions anormales, une methode 
symplectiquc est bien plus cfficaco pour co probleme dc contrSle optimal. L'explication repose sur la theorie de 

I'analysc retrograde en integration numcriquc gcomctriquc. 

Mots-cles : geometrie sous-riemannienne. Martinet, geodesique anormale, integrateur symplectique, analyse 
retrograde 
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Introduction 

This paper is a follow-up to a series of articles that were published in the past decades, see [1, 2j and the 
references therein. There, the authors provide an analytic study of the singularity of the sub-Riemannian 
sphere in the Martinet case. They complement their analysis with numerical computations to represent the 
geodesies and the sphere, and to locate the conjugate points. The integrator used for these computations is an 
explicit Runge-Kutta method of order 5(4). The goal of the present paper is to compare the performances of 
a symplectic integrator versus a non-symplectic one for this optimal control problem. To be more precise, let 
(/7, A, g) be a sub-Riemannian structure where U is an open neighborhood of i?^, A a distribution of constant 
rank 2 and g a Riemannian metric. When A is a contact distribution, there are no abnormal geodesies, and 
a non-symplectic integrator is as efficient as a symplectic one. However, when the distribution is taken as the 
kernel of the Martinet one-form, we show that a symplectic integrator is much more efficient for the computation 
of the normal geodesies and their conjugate points near the abnormal directions. 

Both problems, the Martinet fiat case and a non integrable perturbation, are introduced in Sect. [T] together 
with the corresponding differential equations. Numerical experiments with an expHcit Runge-Kutta method 
and with the symplectic Stormer-Verlet method are presented in Sect. [2] and illustrated with figures. Close to 
abnormal geodesies, the results are quite spectacular. For a relatively large step size, the symplectic integrator 
provides a solution with the correct qualitative behavior and a satisfactory accuracy, while for the same step 
size the non-symplectic integrator gives a completely wrong numerical solution with an incorrect behavior, 
particularly for the non integrable case. The explanation relies on the theory of backward error analysis 
(Sect. [3]). It is related to the geometric structure of the problem and its solutions. 

1 A Martinet type sub-Riemannian structure 

In this section, we briefly recall some results of [1] for a sub-Riemannian structure {U,A,g). Here, U is an 
open neighborhood of the origin in with coordinates q = {x, y, z), and 5 is a Riemannian metric for which a 
graduated normal form, at order 0, is g = (1 -f ay)dx^ + {1 + Px + jy)dy^. The distribution A is generated by 

the two vector fields ^1 = ^ + and F2 — which correspond to A = kerw where uj — dz — ^^dx \s the 
Martinet canonical one-form. To this distribution we associate the affine control system 



where ui{t),U2{t) are measurable bounded functions which act as controls. 

We consider two cases, the Martinet flat case g = dx^ + dy^, an integrable situation, and a one parameter 
perturbation g — dx^ -f (1 + l3x)'^dy'^ for which the set of geodesies is non integrable. 

1.1 Geodesies 

It follows from the Pontryagin maximum principle, see [U [2], that the normal geodesies corresponding to 
g = dx^ -1- (1 + j3xYdy'^ are solutions of an Hamiltonian system 



q ^ Ui{t)Fi{q) + U2{t)F2{q), 



Q = 




(1) 



where q = {x,y, z) is the state, p = {px,Py,Pz) is the adjoint state, and the Hamiltonian is 




In other words, the normal geodesies are solutions of the following equations: 



+ 




Px 




X 



Px 



{l + f3x)3 



y 




Py 




(2) 



Pz 



RR n° 0123456789 



4 



Chyba & Hairer & Vilmart 



Notice that the variables z and Pz do not influence the other equations (except via the initial value pz(0)), so 
that we are actually confronted with a Hamiltonian system in dimension four. For the Martinet flat case (/3 = 0), 
the interesting dynamics takes place in the two-dimensional space of coordinates {y,py). The Hamiltonian is 

where Px and Pz have to be considered as constants. This is a one-degree of freedom mechanical system with a 
quartic potential. Forp^; < < pz, the Hamiltonian H{y,py) has two local minima at {y^±y/—2px/pz, Py = 0), 
which correspond to stationary points of the vector fleld. In this case, the origin (j/ = 0, py = 0) is a saddle point. 

Whereas normal geodesies correspond to oscillating motion, it is shown in [Ull] that the abnormal geodesies 
are the lines z = zq contained in the plane y = 0. For the considered metrics, the abnormal geodesies can be 
obtained as projections of normal geodesies, we say that they are not strictly abnormal. In [2], the authors 
introduce a geometric framework to analyze the singularities of the sphere in the abnormal direction when 
/3 7^ 0. See also [3l|4] for a precise description of the role of the abnormal geodesies in sub-Riemannian geometry 
in the general non-integrable case, i.e., when the abnormal geodesies can be strict. The major result of these 
papers is the proof that the sub-Riemannian sphere is not sub-analytic because of the abnormal geodesies. 

Interesting phenomena arise when the normal geodesies are close to the separatrices connecting the saddle 
point. Therefore, we shall consider in Sect. [2] the computation of normal geodesies with y{0) — and Py{0) > 
but small. 



1.2 Conjugate points 

For the Hamiltonian system ^ we consider the exponential mapping 

which, for flxed qo G , is the projection q{t,qo,po) onto the state space of the solution of starting at 
t = from {qo,Po)- Following the deflnition in [T] we say that the point qi is conjugate to go along q{t) if there 
exists (pojti), ti > 0, such that q(t) = exp^^ ((po) with qi = exp^^ (^(po), and the mapping exp^^ is not an 
immersion at po. We say that qi is the flrst conjugate point if ti is minimal. First conjugate points play a 
major role when studying optimal control problems since it is a well known result that a geodesic is not optimal 
beyond the flrst conjugate point. 

For the numerical computation of the flrst conjugate point, we compute the solution of the Hamiltonian 
system |[T]) together with its variational equation, 

y - J-'VH{y), ^ = J-^V^H{y)-^. (3) 

Here, y = {q,p) and J is the canonical matrix for Hamiltonian systems. It can be shown that for Runge-Kutta 
methods, the derivative of the numerical solution with respect to the initial value, = dyn/dyo, is the result 
of the same numerical integrator applied to the augmented system ([3]), see [F, Lemma VI.4.1]. Here, the matrix 

^ ^ / dq/dqo dq/dpo \ 
V dp/dqo dp/dpo J 

has dimension 6x6. The conjugate points are obtained when dq/dpo becomes singular, i.e., det{dq/dpo) = 0. 

2 Comparison of symplectic and non-symplectic integrators 

For the numerical integration of the Hamiltonian system JT]), where we rewrite ^{q,p) ~ Hq{q,p) and 
^{q,p) — Hp{q,p), we consider two integrators of the same order 2: 

• a non-symplectic, explicit Runge-Kutta discretization, denoted Rk2 (see fSj Sect. II. 1.1]), 

9n+l/2 = qn + ^Hp{qn,Pn) Pn+1/2 = Pn - ^Hq{qn , Pn) 

qn+1 = qn + hHp{qn+i/2,P„+l/2) Pn+1 = Pn ~ hHq{qn+i/2 , Pn+l/2) 
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the symplectic Stormer-Verlet scheme (see e.g. fSJ Sect. VI. 3]), 



Pri+l/2 
9n+l 
Pn+1 



Pn 

q-n 



-Hq{qn,Pn+l/2) 



2 (^p(9"'Pn+i/2) + Hp{q 



(5) 



Pn+1/2 - ■^Hq{qn+l,Pn+l/2) 



where g„ = (x„, ?/„, z„) and p„ = (p 

a:,n i Py,n ■} Pz,n 

). Here, « q{nh), p„ « and h is the step size. 

For the computation of the conjugate points, we apply the numerical methods to the variational equation ([3]). 
Notice that only the partial derivatives with respect to po have to be computed. Conjugate points are then 
detected when det{dqn/dpo) changes sign. We approximate them by linear interpolation which introduces an 
error of size 0{h'^). This is comparable to the accuracy of the chosen integrators which are both of second order. 

Remark 2.1 The Stormer-Verlet scheme |[5]) is implicit in general. A few fixed point iterations yield the 
numerical solution with the desired accuracy. Notice however that the method becomes explicit in the Martinet 
flat case /3 = 0. One simply has to compute the components in a suitable order, for instance Px,n+i, Pz,n+i, 

Py,n+l/2, Un+li Xn+1-, ^n+l, Py,n+1- 

2.1 Martinet flat case 

We consider first the fiat case /3 = in the Hamiltonian system ((2]). As initial values we choose (cf. [T]) 

x(0) = 2/(0) = z(0) = 0, p,(O) = cos0o, P;,(0) ^sin^o, ^^(0) = 10, where 9^ = n - Ur^ (6) 
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Figure 1: Trajectories in the (x,?/)-plane for the fiat case P = 0. 
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so that we start close to an abnormal geodesies, and we integrate the system over the interval [0,9]. 

Figure [T] displays the projection onto the (a;,y)-plane of the numerical solution obtained with different step 
sizes h by the two integrators. The initial value is at the origin, and the final state is indicated by a triangle. The 
circles represent the first conjugate point detected along the numerical solution, while the stars give the position 
of the first conjugate point on the exact solution of the problem. There is an enormous difference between the two 
numerical integrators. The symplectic (Stormer- Verlet) method |(5]) provides a qualitatively correct solution 
already with a large step size h = 0.1, and it gives an excellent approximation for step sizes smaller than 
h = 0.05. On the other hand, the non-symplectic, explicit Runge-Kutta method (0]) gives completely wrong 
results, and step sizes smaller than 10"^ are needed to provide an acceptable solution. An explanation of the 
different behavior of the two integrators will be given in Sect. [3] below. 
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Figure 2: Phase portraits in the (y,pj,)-plane for the fiat case l3 — 0. 



As noticed in Sect.IH the normal geodesies in the fiat case are determined by a one-degree of freedom 
Hamiltonian system in the variables y and Py. We therefore show in Figure [2] the projection onto the {y,py)- 
space of the solutions previously computed with step size h = 0.05. The exact solution starts at (O,sin0o) 
above the saddle point, turns around the positive stationary point, crosses the Py-axis at (0, — sin 6*0), turns 
around the negative stationary point, and then continues periodically. The numerical approximation by the 
non-symplectic method covers more than one and a half periods, whereas the Stormer- Verlet and the exact 
solution cover less than one period for the time interval [0, 9]. Since the conjugate point is not very sensible with 
respect to perturbations in the initial value for Py, the {y,Py) coordinates of the conjugate point obtained by 
the non-symplectic integrator are rather accurate, but the corresponding integration time is completely wrong. 

Table [T] lists the conjugate time obtained with the two integrators using various step sizes. There is a 
significant difference between the two methods. We can see that with the Stormer- Verlet method ^ a step 
size of order h = 10"^ provides a solution with 4 correct digits. A step size a 100 times smaller is needed to get 
the same precision with the non-symplectic method. 



2.2 Non integrable perturbation 

For our next numerical experiment we choose the perturbation parameter f3 — — lO^^ in the differential equation 
([2]). We consider the same initial values and the same integration interval as in Sect. 12. II The exact solution 
is no longer periodic and, due to the fact that /? is chosen negative, its projection onto the {y,py)-space slowly 
spirals inwards around the positive stationary point (see right picture in Figured]). 

Figures [3] and [4] and Table [T] display the numerical results obtained by the two integrators for the differential 
equation ^ with /? = — lO"*. The interpretation of the symbols (triangles, circles, and stars) is the same 
as before. The excellent behavior of the symplectic integrator is even more spectacular than in the fiat case, 
and the pictures obtained for the Stormer- Verlet method agree extremely well with the exact solution. The 



Table 1: Accuracy for the first conjugate time. 
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Figure 3: Trajectories in the (a;, y)-plane for the non integrable case /? = —10 
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Figure 4: Phase portraits in the (?/,Py)-plane for the non integrable case /9 = — 10 



non-symplectic method gives qualitatively wrong solutions for step sizes larger than h — 0.01. In the {y,Py)- 
space it alternatively spirals around the right and left stationary points whereas the exact solution spirals only 
around the positive stationary point. In contrast to the Martinet fiat case, the conjugate point obtained by the 
non-symplectic method is here wrong also in the {y,py)-space. 



2.3 An asymptotic formula on the first conjugate time in the Martinet fiat case 

Now that we have shown the efficiency of symplectic integrators, we can make more precise the asymptotic 
behaviour studied in For the initial values of ^ and /? = 0, consider the ratio 

3X(fc)' 

where ti is the first conjugate time for the normal geodesic, and K{k) is an elliptic integral of the first kind, 

K{k)^ =du, fc = sin(6'o/2). 

"'0 Vl — A;2 gj]^2 y 

By studying analytic solutions for the normal geodesies, it is proved in [l] that this ratio satisfies the inequality 
2/3<i?<l. It follows from a rescaling of the equations iJH) that R is independent oi Pz- 

In Figure m we represent the values of 1 — i? as a function of e = tt — ^o, for various initial values 6*0. The 
numerical results indicate that the ratio R depends on 6*0, and R — > 1^ slowly for 6*0 — > . 
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2 • 10-2 



10-2 



10-'' 10-5 10-4 10-^ 10-2 10-1 10° 
Figure 5: Illustration of the asymptotic behaviour of R (Stormer-Verlet scheme with step size h = IQ-^). 

3 Backward error analysis 

The theory of backward error analysis is fundamental for the study of geometric integrators and it is treated 
in much detail in the monographs of Sanz-Serna & Calvo Hairer, Lubich & Wanner O Chap. IX], and 
Leimkuhler & Reich [^. It allows us to explain the numerical phenomena encountered in the previous section. 

3.1 Backward error analysis and energy conservation 

We briefly present the main ideas of backward error analysis for the study of symplectic integrators, see [U 
Chap. IX]. Consider a system of differential equations 

y^fiy), y{0)^yo (7) 

and a numerical integrator y„+i = ^hiVn) of order p. The idea is to search for a modified differential equation 
written as a formal series in powers of the step size h, 

y = m - fiy) + hFU+i{y) + /i^'+Vp+sly) + • ■ ■ , (8) 

such that yn = y{tn) for tn — nh, n — 0, 1, 2, . . ., in the sense of formal power series. The motivation of this 
approach is that it is often easier to study the modifled equation |[8]) than directly the numerical solution. 

What makes backward error analysis so important for the study of symplectic integrators is the fact that, 
when applied to a Hamiltonian system y — J^^\7H{y), the modified equation ^ has the same structure 
y = J~^'VH{y) with a modified Hamiltonian 

H{y) = H{y) + WHp+,{y) + W+^H^+^^y) + • ■ ■ • 

However, the series usually diverges, so a truncation at a suitable order N{h) is necessary, 

H{y) = H{y) + h^H^+i^y) + ... + h^^-'H^^y). 

This truncation induces an error that can be made exponentially small, by choosing N{h) ^ C/h, see O 
Theorem IX. 8.1]. More precisely, we have that for i„ = nh and /i ^ 0, 

i7(2/„) = i7(yo) + 0(i«e-''°/"). (9) 

as long as the numerical solution {yn} stays in a compact set. On intervals of length 0{e^°/'^^), the modified 
Hamiltonian H(y) is thus exactly conserved up to exponentially small terms. 

3.2 Backward error analysis for the Martinet problem 

Symplectic integrators are successfully appHed in the long-time integration of Hamiltonian systems, for instance 
in astronomy (e.g. the Outer Solar System over 100 million years fH Sect. 1.2.4]), or in molecular dynamics 
O Chap. 11]. Here the situation is quite different because we are interested in the numerical integration of 
Hamiltonian systems on relatively short time intervals. 
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3.2.1 Martinet flat case 

Consider the Martinet problem ^ in the fiat case = 0. Its interesting dynamics takes place in the {y,Py) 
plane, and it is not influenced by the other variables (only by their initial values). We put rj = {y,py), and 
we denote by f{r]) the Hamiltonian vector fleld composed by the corresponding two equations of ([2]). For a 
numerical integrator of order p — 2, the associated modifled differential equation has the form 

?i=f{?i) + h'Mv) + o{h'). (10) 

Consider first the symplectic Stormer-Verlet method. It follows from Sect. 13. II that its modified differential 
equation is Hamiltonian, and from ^ that the modified Hamiltonian H{ri) is preserved up to exponentially 
small terms along the numerical solution. This implies that the numerical solution remains exponentially close 
to a periodic orbit in the {y,py)-space. The critical point {y = 0, Py = 0) is a saddle point also for the modified 
differential equation (because the origin is stationary also for the numerical solution and thus for the modified 
equation). Therefore, any numerical solution starting close to the origin has to come back to it after turning 
around one of the stationary points. The minimal distance to the origin will always stay the same (see the zoom 
in Figure [2]). This explains the good behavior of symplectic integrators. 

For the non-symplectic integrator, the term h^fsir]) is not Hamiltonian. Therefore the solution of the 
modified differential equation (and hence also the numerical solution) is no longer periodic. In fact, it spirals 
outwards and after surrounding the first stationary point, the numerical solution does not approach the saddle 
point sufficiently close, which induces a faster dynamics as can be observed in Figures [T] and [H This causes a 
huge error, because close to the saddle point the numerical solution is most sensible to errors. 

3.2.2 Non integrable perturbation 

In this case, the argument in the comparison of symplectic and non-symplectic integrators is very similar to 
the discussion of the Van der Pol's equation in [SJ Sect. XII. 1]. For /? 7^ (non integrable perturbation), the 
dynamics takes place in the four dimensional space with variables t] = {x,y,Px,Py)- In this space the system 
^ becomes 

v = fiv) + I3g{v) 

where f{r]) is the Hamiltonian vector field corresponding to /? = and g{y) = C(l) depends smoothly on (3. 
Here, the modified equation becomes 

v = f{v) + Pgiv) + h'fsiv) + Oie + ph"), 

where the perturbation term h? fzijf) is the same as for the Martinet flat case. 

For the symplectic integrator, the perturbation l3g{r]) has the same effect for the original problem as for 
V = f iv) + h'^ hiv) + ■ ■ ■ ■ This explains the correct qualitative behavior for small h and small f3. There is no 
restriction on the step size h compared to the size of /?. 

For the non-symplectic integrator, each of the perturbation terms l3g(r]) and h? fziif) destroys the periodic 
orbits in the subsystem for the {y,Py) variables, and the dominant one will determine the behavior of the 
numerical solution. Only when /i^ <C \(3\, the numerical solution will catch the correct dynamics of the problem. 
In Figures [3] and m where (3 = —10"'*, this condition is not satisfied for h > 10~^. Since /? is chosen small 
and negative, the two perturbation terms are conflicting. The term l3g{ri) causes the solution to spiral around 
the positive stationary point, whereas the term h? fzijf) causes it to spiral alternatively around both stationary 
points. For too large step sizes the qualitative behavior of the non-symplectic integrator ^ is thus completely 
wrong. 

Remark 3.1 The problem ^ with (3 — has a lot of symmetries. In the (?/,Pj^)-space the orbits are symmetric 
with respect to the y-axis and also with respect to the Pj,-axis. If we apply a symmetric numerical integrator (not 
necessarily symplectic) , it is possible to prove the same qualitative behavior as for the symplectic Stormer-Verlet 
method. This follows from the fact that the solution of the modified equation (numerical orbit) corresponding 
to a symmetric method has the same symmetry properties as the exact fiow (see [U Sect. IX. 2] for precise 
statements). Consequently, in the {y,Py) plane and for 13 — 0, the solution will stay exponentially near to 
a closed orbit, as it is the case for symplectic integrators. In the non integrable case, the good behavior of 
symmetric methods can be explained as in Sect. 13.2.21 for symplectic methods. 
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